function constrains = constrains_C0(dof_map, TE, ET, FE_Order)
% function constrains = constrains_C0(dof_map, TE, ET, FE_Order)
% get the global smoothness conditionsidx_mat = zeros(3,d,3); 
% the column idx 1 2 3 -1 -2 -3
%
% The returned value is orginized as in
%                      
%
idx_c0 = zeros(FE_Order + 1, 3);
cr = cr_pattern(FE_Order);
for i = 1:3
    idx_c0(:,i) = cr_indices(0, FE_Order, i, cr);
end
E_in = find(ET(:,2)~=0); nE_in = length(E_in);
n_row = nE_in*(FE_Order + 1);
constrains = zeros(n_row, 2);
pos = 0;
for k = 1:nE_in
    eg = E_in(k);
    tri = ET(eg,1); nei = ET(eg,2); % find two triangles who have eg.
    eg2tri = find(TE(tri,:)==eg);  % local index of eg in tri
    eg2nei = find(TE(nei,:)==eg); % local index of nei in tri
    row_idx = pos + (1:(FE_Order + 1));
    % first set c0 condition:
    
    constrains(row_idx, 1) = dof_map(idx_c0(:,eg2tri), tri);
    % set the row index for nei:
    if mod(eg2tri + eg2nei,2)~=1  % the inverse direction when any of them are 2, good criteria!!!
        row_idx = pos + ((FE_Order + 1):-1:1);
    end
    constrains(row_idx, 2)  = dof_map(idx_c0(:,eg2nei), nei);
    
    pos = pos + FE_Order + 1;
end

% % need to get rid of boundary constrain, which have been considered as
% % dirichlet boundary conditions
% flag = zeros(n_row, 1);
% [row,col] = size(dof_boundary);
% for irow = 1: row
%     for jcol = 1: col
%         flag = flag + (constrains(:,1) == dof_boundary(irow, jcol)) + (constrains(:,2) == dof_boundary(irow, jcol));
%     end
% end
% % find all the necessary boundary constrains
% idx = find(flag ~= 0);   % !!!!  THIS CONDITION IS NOT SUFFICIENT !!!
% % and obtain them purely
% constrains_new = constrains(idx, :); 
% n_row = length(idx);

% H = sparse((1:n_row)'*[1, 1], constrains, [ones(n_row,1),  -1*ones(n_row,1)], n_row, max(max(dof_map)));

end

